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We introduce a one-dimensional model for a CuO chain, with holes and S = 1/2 spins localized in 
3da,2_^2 orbitals, and pa oxygen orbitals without holes in the ground state. We consider a single 
hole doped at an oxygen site and study its propagation by spin-flip processes. We develop the 
Green’s function method and treat the hole-spin coupling in the self-consistent Born approximation, 
similar to that successfully used to study polarons in the regular t-J model. We present an analytical 
solution of the problem and investigate whether the numerical integration is a good approximation 
to this solution. 
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Soon after the discovery of high-Tc superconductors 
(HTS), it was realized that the t-J model could be re¬ 
garded as the prime candidate for the theoretical de¬ 
scription of the phenomenon [1]. A deeper insight into 
the dynamics of a hole propagation in a Mott insulator 
with antiferromagnetic (AF) order might yield a clue 
to the understanding of the origin of superconductiv¬ 
ity in cuprates [2]. The effective t-J model arises from 
mapping the three-band p-d model [3, 4] onto the cop¬ 
per dx 2 _y 2 states and provides the simplest description 
of the electronic states in Cu02 planes of HTS. Recent 
studies indicate, however, that oxygen p-orbital states 
strongly renormalize the quasiparticle (QP) energy, both 
for antiferromagnetic [7] and ferromagnetic [8] systems. 
Importance of oxygen orbitals may be recognized both 
from the prominent role played in cuprates by the second 
neighbor effective hopping P, for instance in the stability 
of stripe phases [5] , and from numerical studies comparing 
hole-doped and electron-doped systems [6]. 

These circumstances have led us to study the polaron 
dynamics in the extended t-J model which includes the 
oxygen states. In analogy to the well known Cu02 sys¬ 
tem, initially the 2p oxygen states are completely filled 
with electrons, while the < 7 ^, 2_^2 states are half-filled and 
hosting localized S = 1/2 spins. The superexchange J 
couples the spins on Cu sites [9] and stabilizes AF order 
in the ground state. In the charge-transfer insulator there 
are no itinerant charges in the p states. Therefore, to 
activate any kinetic energy or interaction whatsoever, one 
needs to inject charges into the p states, which is usually 
achieved by means of hole doping La 2 -a;(Sr,Ca)a,Cu 04 , or 
by substitutional transition metal impurities of different 
valence than Cu ions. Here we consider the case of a 
single hole injected into the p-band. 

The complexity of hopping processes in a realistic 
CuOs chain in YBa 2 Cu 307 [10] motivates a simpler one¬ 
dimensional (ID) model of a CuO chain, see Fig. 1. The 
advantage of this linear structure is that interoxygen hop¬ 



FIG. 1. Schematic representation of the CuO model with 
the conventions used throughout this paper. The dashed line 
frame indicates the magnetic unit cell in real space. 


ping is avoided and one may focus on the hole dynamics 
of a carrier injected into one of the oxygen orbitals, which 
follows from three-site processes, as shown below. Such 
processes are of crucial importance in systems where other 
hopping processes cannot occur, as in t 2 g systems with 
Ising-like superexchange, where they cause weak hole 
propagation [11] and stabilize bond stripes [12], and in 
the quantum compass model [13]. 

We formulate the simplest ID t-J-like model with hole 
dynamics in p orbitals. The model Hamiltonian, 


H = HoPV, ( 1 ) 

consists of magnetic Heisenberg exchange (2a), enforcing 
the AF ground state in the “Cu” d states and a three-site 
p-d-p spin-flip hopping (2b), obtained in second order 
perturbation theory, acting as the interaction: 

no = jJ 2 {Si-Si+ 2 i-S^), (2a) 

i 

V = t + H-C.), (2b) 

i,a 

where ^ = a/4, a is the magnetic unit cell length, the 
constants J (superexchange) and t (hopping) are defined 
as positive, and the summation extends over all the crystal 
unit cells. The model follows in perturbation theory from 
the respective itinerant model for the CuO chain. Note 
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that this simple Hamiltonian does not include any free 
electron kinetic terms which are blocked due to large 
Coulomb repulsion on dx^_y 2 orbitals, and therefore the 
free electrons are dispersionless. For a graphic illustration 
of the model and the various conventions pertaining to it, 
consult Fig. 1. 

To solve this toy model in the self-consistent Born 
approximation (SCBA) [14], one should first notice that 
the pairs of fermion operators at Cu sites are equivalent 
to spin-flip processes. 


d,,dl = -Sf, (3) 

where the right hand side is the spin raising/lowering op¬ 
erator, depending on the value of the spin index a. SCBA 
is typically employed using the Linear Spin Wave (LSW) 
theory, where one expands the spin operators in powers 
of bosonic operators up to the second order using the 
Holstein-Primakoff representation. All the calculations 
are done here in the magnetic unit cell, as the expansion 
of (2b) around the AF ground state is irreducible to the 
crystal unit cell. 

After performing the standard steps of Fourier and 
Bogoliubov transformations in the two sublattice frame¬ 
work, one arrives at the following representation of the 
Hamiltonian (2): 


Uo = J2uig{alag + /3l/3g), 


(4a) 


^ ^ X] ^'=9 + /5?) + H-c. , (4b) 
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expression for the first order self-energy: 
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where I 2 is a 2 x 2 identity matrix to account for the 
aforementioned spin degrees of freedom and the summa¬ 
tion extends over the first Brillouin Zone (BZ). The above 
equation can be solved analytically for the present ID 
case, yielding 
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however the details of this derivation are beyond the scope 
of this article. 

On the other hand, the self-energy (8) can also be 
calculated numerically by simple lattice summation over 
the BZ. After this is accomplished, E^^^(/c,co’) can be used 
to obtain the first order Green’s function via the equation 


G{k,uj) = [{uj \ (10) 


which forms the basis of the self-consistent solution. Equa¬ 
tion (10) can then be inserted into (6) to obtain the second 
order self-energy, and so on. However, usually the first 
order solution suffices and further iterations do not change 
the result significantly. Once the Green’s function is cal¬ 
culated, one can extract the physical information in the 
form of the spectral function 

^(fc,w) =-/-9[TrG'(fc,w)], (11) 


where ujq = J|sin(g'/2)| is the magnon dispersion relation 
in the ID AF model (2a), N is the number of magnetic 
unit cells in the system, and 

n, = sin^sin|(ug + t)g) (5) 

is the electron-magnon vertex function. The opera¬ 
tors {aq,/3q} are the bosonic Bogoliubov operators and 
{uq,Vq} are the Bogoliubov transformation coefficients. 
More details on the above derivation can be found in the 
literature concerning polarons in the t-J model, particu¬ 
larly in [14]. 

To obtain the SCBA Green’s function solution, one 
needs to calculate the self-energy, which then serves as 
the basis for the first order Green’s function: 

E{k,uj) = {0\p,,vg{uj)Vpl\0), (6) 

Q{uj) = [uj ^if] - H]“\ (7) 

where |0) is the AF ground state and Q(uj) is the resolvent, 
or the Green’s function operator. Note that, because of 
the spin degrees of freedom, 'E(k,uj) is in fact a 2 x 
2 matrix, although only in a trivial, diagonal manner. 
Inserting (4) into (6) one quickly arrives at the following 


which has a direct relation to photoelectron spectroscopy 
experiments. 

Here we will present the numerical results (for self¬ 
energy in the first order) obtained by the method de¬ 
scribed above, and contrast them against the first order 
analytic result, as mentioned earlier. As can be seen 
from (8), the natural unit of energy for this problem is 
the parameter J, therefore the results are plotted using 
the energy uj/J. 

Figure 2 shows the density map of the numerically 
obtained spectral function A{k^uj) (in gray-scale) for the 
problem under consideration with t = 1. Bear in mind 
that in order to expose the low-amplitude features of 
the spectrum, a nonlinear tanh-scale has been employed, 
hence the very strong “broadening” of the bands. 

On top of the numerical result we have plotted the 
analytic solution (solid red line). This plot was obtained 
by solving the equation 

uj — 5R[E^^^(/c,cl; + irj)] = 0, (12) 

which defines the location of the QP bands, sans the 
broadening of the spectrum, i.e., the resulting curve cor¬ 
responds to the QP maxima of (11). Next we will discuss 
the differences between the two results. 
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FIG. 2. Numerically calculated spectral function A{k,(jj) (gray 
scale) and the analytic solution (solid red lines) for t = 1 and 
rj = 0.01. Note the nonlinear tanh-scale. 

First of all, one notices that both the numerical and 
the analytic results display two solution branches: (i) the 
lower branch with negative energy which is a bound QP 
state (since the reference free particle dispersion is 0), and 
(a) the upper branch, which is an excited QP state — in 
the present case, as the analytic result seems to suggest, 
it exists only for k ^ 7r/4. 

It is evident from Fig. 2 that in the case of the lower 
branch the analytic result reproduces the numerical so¬ 
lution quite well in the whole range of k values. On the 
other hand, one can easily notice a sharp discrepancy 
between the two results in the case of the upper branch. 
While the analytic result ends rather abruptly around the 
energy cj = 1, the numerical solution extends well below 
this point. Since the analytic solution does not employ 
any approximations (beyond the physical ones), we rather 
expect the numerical approach to be unsound. 

It is our supposition that the problem lies within the 
too naive an approach to the numerical integration. More 
specifically, the function being integrated in (8) is propor¬ 
tional to {uj — |sin(g'/2)|)“^, which has a singularity line 
in the range uj G [0,1]. Integrating such a function numer¬ 
ically using simple quadrature rules, such as the rectangle 


rule in this case, introduces noticeable systematic errors 
into the solution. Such errors should be especially visible 
in the range in which the divergence occurs. As can be 
seen in Fig. 2, this is indeed the case since the discrepancy 
between the two solutions is most pronounced exactly in 
the range where the singularities occur. 

Further, from (8) one can see that in order to calculate 
S(/c,ci;), one needs to integrate all the (k — q,uj — ujq) 
elements, which means that the errors introduced in the 
lower energies propagate to the solution at higher energies, 
which would explain the mismatch between the solutions 
also for cj > 1. Moreover, this implies that a naive 
numerical solution of the problem presented herein is 
reliable only until the energy cj 0. 

In summary, we have presented the exact analytic result 
and an approximate numerical solution of the formulated 
problem. By comparing them we have emphasized that 
an unsophisticated approach to numerical methods can 
be detrimental to the reliability of the SCBA results. 
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